% PCA main program
% Analizes an image sequence, performs PCA
% EGC 2009/04/15

%% 1st sequence (Makes a movie of fluorescent images)
addpath(genpath('D:\Users\Turtoise\'));
path1 = uigetdir('D:\Users\Data\IVIS\Imaging Session 7 04-22-2009\AEK20090422155922_SEQ\','Sequence folder');
if path1
    [X,k,t,D,F,R,I,photo] = make_movie(path1,'nodisp',0.1);
%     [X,k,t] = make_movie(path1);
else
    disp('Cancelled by user');
end

%% Example of image normalization
R_temp = R{randi(length(R),1)};         % A random image
% Background and flatfield correction
I_temp = double(R_temp - D) ./ double(F - D);
M = F - D;
M = mean(M(:));                         % Determine the average pixel value M in the corrected flat frame (F 	 D).
I_temp = M * I_temp;                    % Corrected Image

figure
colormap gray
set (gcf,'color','white');

subplot(2,2,1)
imagesc(R_temp)
axis square
title('Raw data','FontSize',12)

subplot(2,2,2)
imagesc(F)
axis square
title('Fluorescent Reference','FontSize',12)

subplot(2,2,3)
imagesc(D)
axis square
title('Bias image','FontSize',12)

subplot(2,2,4)
imagesc(I_temp)
axis square
title('Corrected image','FontSize',12)

clear R_temp I_temp

%% Example of Mean Fluorescence Overlay
threshold = 0.2;

h1 = figure;
set(h1,'visible','off')                         % Invisible window
set (h1,'color','white');
moyen = reshape(mean(X,2),[256 256]);
% moyen = imadjust(moyen,stretchlim(moyen),[]);     % Don't use
moyen = moyen ./ max(moyen(:));
moyen = moyen .* (moyen>threshold);
imagesc(moyen)
colormap hot
hax1 = gca;
cbar_axes = colorbar;
axis square
axis off
title('Mean Fluorescence','FontSize',14)
carte_orig = get(gcf,'Colormap');
carte = ind2rgb(gray2ind(moyen),carte_orig);
close(h1)                                       % Close invisble window

h2 = figure;
set (h2,'color','white');
hax2 = imagesc(photo);
colormap gray
axis square
hold on
hax3 = imshow(carte);
hax4 = gca;
title('Mean Fluorescence Overlay','FontSize',14)
set(hax3, 'AlphaData',moyen>threshold);

clear h1 h2 hax1 hax2 hax3 hax4 threshold carte carte_orig cbar_axes

%% Principal Component Analysis
% PCA on original data princomp(X)
% [PC,SCORE,latent,tsquare] = princomp(X);
% PCA on centered and scaled data princomp(zscore(X));
[PC,SCORE,latent,tsquare] = princomp(zscore(X));
r = reshape(SCORE(:,2),[256 256]);
g = reshape(SCORE(:,3),[256 256]);
b = reshape(SCORE(:,4),[256 256]);

rneg = -r.*(r<0);
rpos = r.*(r>=0);
gneg = -g.*(g<0);
gpos = g.*(g>=0);
bneg = -b.*(b<0);
bpos = b.*(b>=0);

% -------------------------------------------------------------------------
% Normalize between 0 and 1 (optional, to use instead of imadjust)
% r = r - min(r(:));
% r = r./max(r(:));
% g = g - min(g(:));
% g = g./max(g(:));
% b = b - min(b(:));
% b = b./max(b(:));
% -------------------------------------------------------------------------

% Eliminates 2% outliers and normalizes values between 0 and 1
r = imadjust(r,stretchlim(r),[]);
g = imadjust(g,stretchlim(g),[]);
b = imadjust(b,stretchlim(b),[]);

% Eliminates 2% outliers and normalizes values between 0 and 1
rneg = imadjust(rneg,stretchlim(rneg),[]);
gneg = imadjust(gneg,stretchlim(gneg),[]);
bneg = imadjust(bneg,stretchlim(bneg),[]);

% Eliminates 2% outliers and normalizes values between 0 and 1
rpos = imadjust(rpos,stretchlim(rpos),[]);
gpos = imadjust(gpos,stretchlim(gpos),[]);
bpos = imadjust(bpos,stretchlim(bpos),[]);

% Merging components
% Principal components merged image
pci = reshape([r g b],[256 256 3]);
figure; image(pci); set (gcf,'color','white'); axis square;
title('RGB merged image (PC2,PC3,PC4)','FontSize',18)

% Negative Principal components merged image
pcineg = reshape([rneg gneg bneg],[256 256 3]);
figure; image(pcineg); set (gcf,'color','white'); axis square;
title('Negative RGB merged image (PC2,PC3,PC4)','FontSize',18)

% Positive Principal components merged image
pcipos = reshape([rpos gpos bpos],[256 256 3]);
figure; image(pcipos); set (gcf,'color','white'); axis square;
title('Positive RGB merged image (PC2,PC3,PC4)','FontSize',18)

%% Plot PCA time courses
figure
set (gcf,'color','white');
h = plot(t,PC(:,1),'-k.',t,PC(:,2),'-r.',t,PC(:,3),'-g.',t,PC(:,4),'-b.');
xlim([t(1) t(end)]);
set(gca,'FontSize',16)
set(h,'LineWidth',2,'MarkerSize',16)
    legend('PC1','PC2','PC3','PC4','Location','NorthWest');
xlabel('Time (s)','FontSize',18)
title('PCA time courses','FontSize',18)

% Accumulated variance
% v = latent ./ sum(latent);
% for k=1:length(latent)
%     v2(k) = sum(v(1:k));
% end

%% Transparency Mapping

% photo = imread([path2 '\photograph.TIF']); % Photo Image
% photo = imresize(photo,[256 256]);
figure
set (gcf,'color','white');
% imagesc(photo)
% colormap gray
% axis square; axis off

subplot(2,2,1)
imagesc(photo)
colormap gray
axis square
bg = zeros([size(photo) 3]);
PC1 = reshape(SCORE(:,1),[256 256]);
PC1 = imadjust(PC1,stretchlim(PC1),[]);
bg(:,:,1) = PC1;
% bg(:,:,2) = PC1;
bg(:,:,3) = PC1;
hold on
PC1 = reshape(SCORE(:,1),[256 256]);
PC1 = imadjust(PC1,stretchlim(PC1),[]);
hm = imshow(bg);                   % Coefficients correspond to mean image
title('PC1 (corresponds to the mean image)','FontSize',12)
hold off

subplot(2,2,2)
imagesc(photo)
colormap gray
axis square
red = zeros([size(photo) 3]);
red(:,:,1) = r;
hold on
hr = imshow(red);
title('PC2','FontSize',12)
hold off

subplot(2,2,3)
imagesc(photo)
colormap gray
axis square
green = zeros([size(photo) 3]);
green(:,:,2) = g;
hold on
hg = imshow(green);
title('PC3','FontSize',12)
hold off

subplot(2,2,4)
imagesc(photo)
colormap gray
axis square
blue = zeros([size(photo) 3]);
blue(:,:,3) = b;
hold on
hb = imshow(blue);
title('PC4','FontSize',12)
hold off

% Composite image
% figure
% h = imshow(pci);

% Use PC as the AlphaData for the solid PCA colors, overlaying on photo
% set(h, 'AlphaData', reshape(SCORE(:,1),[256 256]) ./ max(SCORE(:,1)))
% threshold = 0.99;
threshold = min([min(r(:)) min(g(:)) min(b(:))]);
set(hr, 'AlphaData',pci(:,:,1).*(pci(:,:,1)>threshold));
set(hg, 'AlphaData',pci(:,:,2).*(pci(:,:,2)>threshold));
set(hb, 'AlphaData',pci(:,:,3).*(pci(:,:,3)>threshold));
set(hm, 'AlphaData',PC1.*(PC1>threshold));

% -------------------------------------------------------------------------
%% Create mask from PC1
mask = PC1>0;

%% Create rectangular ROI (for each region you want to eliminate)
figure
h_im = imagesc(mask); axis square; set (gcf,'color','white'); colormap gray
e = imrect;

%% Execute once for each ROI removal
BW = createMask(e,h_im);
figure; imagesc(BW); axis square; set (gcf,'color','white'); colormap gray
mask = bitor(mask,BW);

%% Show final mask
delete(e);                                  % Deletes ROI
% close all
figure
h_im = imagesc(~mask);                      % Show negative for better visualization
axis square; set (gcf,'color','white'); colormap gray
title('Mask','FontSize',18)

%% Apply mask to all sequence
nmask = repmat(mask(:),[1 k]);
X = X.*(~nmask);

%% Run PCA only on the mouse new data
% Execute from "Principal Component Analysis" cell
% to "Use PC as the AlphaData for the solid PCA colors" cell
